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ABSTRACT 

We study the soft excess variability of the narrow line Seyfert 1 galaxy IRAS 13224-3809. We considered all five archival XMM- 
Newton observations, and we applied the ‘flux-flux plot' (FFP) method. We found that the flux-flux plots were highly affected by 
the choice of the light curves’ time bin size, most probably because of the fast and large amplitude variations, and the intrinsic non¬ 
linear flux-flux relations in this source. Therefore, we recommend that the smallest bin-size should be used in such cases. Hence, we 
constructed FFPs in 11 energy bands below 1.7 keV, and we considered the 1.7-3 keV band, as being representative of the primary 
emission. The FFPs are reasonably well fitted by a ‘power-law plus a constant’ model. We detected significant positive constants 
in three out of five observations. The best-fit slopes are flatter than unity at energies below ~ 0.9 keV, where the soft excess is 
strongest. This suggests the presence of intrinsic spectral variability. A power-law-like primary component, which is variable in 
flux and spectral slope (as T oc AjJ/) and a soft-excess component, which varies with the primary continuum (as F exC ess E 94 ^,), 
can broadly explain the FFPs. In fact, this can create positive ‘constants’, even when a stable spectral component does not exist. 
Nevertheless, the possibility of a stable, soft-band constant component cannot be ruled out, but its contribution to the observed 0.2-1 
keV band flux should be less than ~ 15%. The model constants in the FFPs were consistent with zero in one observation, and negative 
at energies below 1 keV in another. It is hard to explain these results in the context of any spectral variability scenario, but they may 
signify the presence of a variable, warm absorber in the source. 
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1. Introduction 


It is generally accepted that AGN are powered by accretion of 
matter onto a central supermassive black hole (SMBH) of mass 
Mbh ~ 10 fi 9 M 0 . The matter is thought to accrete in a disc that 
is geometrically thin and optically thick (Shakura & Sunyaev 
1973). If the released gravitational energy is dissipated locally 
on the disc, then we expect a multi-temperature, blackbody emis¬ 
sion component to dominate the spectral energy distribution of 
these objects, with a maximum temperature of ~ 10 5 K that peaks 
in the ultraviolet. Indeed, observations of AGN suggest the pres¬ 
ence of a peak in the optical-UV continuum known as the ‘big 
blue bump’ (Malkan & Sargent 1982; Laor et al. 1997; Zheng 
st al. 1997; Shang et al. 2005). 

AGN are strong X-ray emitters. These X-rays are believed 
to be triggered by Compton up-scattering of the disc photons off 
hot electrons surrounding the disc, in a hot (~ 10 9 K) medium, 
usually referred to as the X-ray corona. The X-ray photons are 
emitted isotropically, and although the geometry of disc-corona 
is currently unknown in AGN, they may reflect off the accre¬ 
tion disc and/or the dusty torus (0.1-10 pc from the SMBH) , 
thereby producing an additional ‘reflection’ component in the 
X-ray spectrum (e.g. George & Fabian 1991). 

AGN X-ray spectra at energies above ~ 2 keV have a power- 
law-like shape. The extrapolation of the best-fit power-law mod¬ 
els to energies lower than 2 keV reveals in many AGNs a spectral 
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component in excess of the extrapolated hard-band continuum. 
This excess emission is know as the soft X-ray excess. This com¬ 
ponent was discovered over 30 years ago (Singh et al. 1985; Ar- 
naud et al. 1985), and its source has been debated ever since. 

Originally, it was suggested that soft X-ray excess represents 
the high-frequency tail of the disc emission in AGN. In fact, the 
soft excess in many AGN could be fitted by a blackbody model 
with a best-fit temperature in the range 0.1 - 0.2 keV (Walter 
& Fink 1993; Czerny et al. 2003), however, this temperature is 
significantly higher than the maximum temperatures expected in 
AGN accretion discs. It has also been proposed that the soft ex¬ 
cess could arise from Compton up-scattering of disc photons, 
in a ‘warm’ medium of an electron population with a temper¬ 
ature much lower and an optical depth much higher than those 
of the X-ray corona that are responsible for the emission at en¬ 
ergies above 2keV (e.g. Czerny & Elvis 1987; Magdziarz et al. 
1998; Janiuk et al. 2001). However, Gierlinski & Done (2004) 
show that the temperature associated with this Comptonisation 
region is constant over a wide range of AGN luminosity and 
black hole mass, a result that requires a fine-tuning between the 
optical depth and the heating/cooling ratio of this region. They 
later argue (Gierlinski & Done 2006) that a “smeared” absorp¬ 
tion model can provide an explanation for this component, but 
subsequent studies show that line-driven AGN accretion disc 
winds cannot reproduce the soft excess (Schurch & Done 2007, 
2008). 

Recently, Done et al. (2012) have revived the original idea of 
the soft excess arising from the disc itself. They propose a model 
where the gravitational energy released by accretion at small 
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radii is split between powering optically thick, Comptonized disc 
emission, which is responsible for the soft excess, and an opti¬ 
cally thin corona above the disc, which is responsible for the 
high-energy X-ray emission. Finally, X-ray reflection of the in¬ 
ner disc by the X-ray source can also account for the soft excess. 
If the accretion disc is mildly ionized and the X-rays illuminate 
its innermost region, excess emission below ~ 2 keV could ap¬ 
pear because of line and bremsstrahlung emission from the hot 
disc layers. Crummy et al. (2006) explain the soft excess in many 
AGN using a relativistically-blurred version of the disc reflection 
model of Ross & Fabian (2005), but recent studies have shown 
that this model may not be sufficient in some cases (see Lohfink 
et al. 2012). 

In this work, we apply the flux-flux plot (FFP) method using 
archival XMM-Newton data of the narrow line Seyfert 1 (NLS1) 
galaxy IRAS 13224-3809. Our aim is to study the soft X-ray 
variability in this source, and constrain the origin of its soft X- 
ray excess in a model-independent fashion. The FFP method has 
recently been applied by Noda et al. (2011, 2013) to a hand¬ 
ful of X-ray bright quasars with the same goals. They detect a 
soft component in the 0.5-3 keV band of a few AGN, which is 
less variable than the primary source, and they interpreted this 
in terms of thermal Comptonization of the disc emission. This 
method was first developed by Churazov et al. (2001) to iden¬ 
tify variable and stable components in the X-ray spectrum of 
the black hole binary Cygnus X-l. It was later applied by Tay¬ 
lor et al. (2003) to study the X-ray spectral variability of X-ray 
bright Seyferts and to many other cases since. 

The FFP method is straightforward in its implementation and 
is relatively efficient in detecting spectral components that are 
less variable than the X-ray primary emission. If present, these 
components result in positive constants in the FFP plots, which, 
in general, show a strong positive correlation between the flux in 
various energy bands (for AGN). The FFP variant of the method 
is particularly effective when studying the fast spectral variabil¬ 
ity in AGN on timescales when the flux in various energy bands 
cannot be accurately determined. 

IRAS 13224 - 3809 (z = 0.066) was first detected in X-rays 
during the ROSAT all-sky survey in 1992 as a high-luminous X- 
ray source with an X-ray luminosity Lx - 3 x 10 44 erg/s (Boiler 
et al. 1993). IRAS 13224 - 3809 is highly variable "in X-rays. 
Its spectrum shows a soft X-ray excess below ~ 1.5 ke V (Boiler 
et al. 1996,2003) and strong relativistic effects (Ponti et al. 2010; 
Fabian et al. 2013). Chiang et al. (2015) recently presented a de¬ 
tailed study of the X-ray spectral variability of the source. They 
model its X-ray spectrum in the context of relativistic disc reflec¬ 
tion models. They found that the reflected emission is much less 
variable than the X-ray continuum, and as the source flux drops, 
the spectrum becomes progressively more reflection-dominated. 
The X-ray reflection interpretation is broadly consistent with the 
results of Emmanoulopoulos et al. (2014), who studied the soft- 
band time lags in this source. They found that these are consis¬ 
tent with the hypothesis of a point-like source which is located 
at a height of ~ 3 gravitational radii (r g ) above a ~ 10 7 M G black 
hole (BH). Both the spectral and the timing studies mentioned 
above suggest that most of the soft X-ray emission in this source 
is due to X-ray illumination of the inner disc. 

In this paper, we use all the archival XMM-Newton data of 
IRAS 13224-3809. The fact that it is quite bright and highly 
variable in X-rays, as well as the strong relativistic effects seen in 
this source, make it an interesting target to apply the FFP method 
and study its soft X-ray variability. In Section 2 we present the 
observations and the data reduction, and in Section 3 we present 
the results from the model-fitting of the FFPs. Finally, in the 


Table 1. Log of observations of IRAS 13224 - 3809. The last column 
lists the 1.7-3 keV band flux, assuming a power-law model, taking the 
Galactic absorption into consideration. 


Obs./date 

Net exp. 
(ks) 

Pn-mode 

1.7-3 keV flux 
(10~ 13 erg/s/cm 2 ) 

1/2002-01-19 

61 

FW 

i 94+o.ti 

-0.16 

2/2011-07-19 

120 

FW 

2 43+0.12 
-0.19 

3/2011-07-21 

120 

LW 

1 76 +on 

A *' -0.18 

4/2011-07-25 

110 

LW 

0 94 +0 - 10 
-0.14 

5/2011-07-29 

120 

LW 

1 QC+0.10 

1 * 7 -0.12 


last section we provide a brief summary of our main results and 
discuss their physical implications. 

2. Observations and data reduction 

Five observations of IRAS 13224-3809 by the XMM-Newton 
satellite (Jansen et al. 2001) are available in the XMM-Newton 
Science Archive (XSA). The observation log is listed in Table 1. 
The first observation was done in 2002 January 19 (Obs. ID 
0110890101, hereafter Obs. 1), and the last four observations 
from 2011 July 19 to 29 (Obs. IDs 0673580101, 0673580201, 
0673580301 and 0673580401, hereafter, Obs. 2, Obs. 3, Obs. 4, 
and Obs. 5, respectively). All instruments were working success¬ 
fully during the five observations. We have considered only the 
data provided by the EPIC-pn camera (Strtider et al. 2001). 

The first two observations were taken in full window (FW) 
imaging mode and the others in large window (LW) imaging 
mode. We reduced all data using the XMM-Newton Science 
Analysis System (SAS v. 13.0.0) and the latest calibration files. 
The data were cleaned for high background flares and were se¬ 
lected using the condition PATTERN < 4. The net exposure for 
each observation is also listed in Table 1. Source spectra and light 
curves were extracted from a circle of radius 40". The respec¬ 
tive background spectra and light curves were extracted from an 
off-source circular region located on the same CCD chip, with 
a radius approximately twice that of the source to ensure a high 
signal-to-noise ratio. Pileup was checked and found to be neg¬ 
ligible in all observations. Background-subtracted light curves 
were produced using the SAS task EPICLCCORRfor different en¬ 
ergy bands and time bin sizes. The choice of the energy bands 
and the time binning is discussed in detail in the next section. 

Response matrices were produced using the SAS tasks 
RMFGEN and ARFGEN. Spectral model-fitting was performed with 
XSPEC v. 12.8.1 (Arnaud 1996). We report the lcr errors on the 
model parameter and flux estimates, and 3cr upper/lower limits 
when relevant. 

3. Flux-Flux analysis 

3.1. The ‘continuum’ band 

Our first task was to identify the energy band that is the most 
representative of the X-ray continuum emission. To this end, we 
fitted an absorbed power-law (PL) model to the 3 - 10 keV band 
spectra for all observations (i.e. wabs*powerlaw in XSPEC ter¬ 
minology). We considered only Galactic absorption, which for 
the line of sight towards IRAS 13224 - 3809 corresponds to an 
equivalent hydrogen column density of Nh = 5.34 x 10 2,) cm 2 
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Energy (keV) 

Fig. 1 . Ratio of the observed spectra over a power-law model fitted to 
the 3 - lOkeV band and extrapolated to lower energies (Obs. 1, 2, 3, 4, 
and 5 spectra are shown as black, red, cyan, blue, and magenta points, 
respectively). The two vertical black lines indicate the energy band we 
chose as representative of the X-ray continuum. 


(Kalberla et al. 2005). Figure 1 shows the ratio of the observed 
spectra over the absorbed PL model for the full energy band 
(0.2-10.0 keV) in all observations. Based on the resulting^ 2 val¬ 
ues, a PL cannot be accepted as a good fit model in the 3 -lOkeV 
band. Nevertheless, this model provides a reasonably good base¬ 
line model to account for the hard X-ray emission in the source. 
The extrapolation of the best-fit PL to energies below ~ 1 keV 
clearly reveals a strong soft excess in all observations. Notably, 
the strength of this component is also variable. 

Based on the plot shown in Fig. 1, we chose the 1.7-3.0keV 
band as the most representative of the X-ray ‘continuum’ (indi¬ 
cated by the vertical lines in Fig. 1). The choice of the lower en¬ 
ergy limit was dictated by the requirement to be as low as possi¬ 
ble (to increase the signal-to-noise in the continuum band) and at 
the same time to be the least affected by the soft excess emission. 
The choice of high energy limit was dictated by the possibility 
that the X-ray reflection component may strongly affect the 3- 
lOkeV band in this source (Fabian et al. 2013). In any case, the 
addition of the 3-10keV counts in the continuum band does not 
increase the count rate significantly in the band we chose, while 
it does increase its error rate considerably (because the source 
is relatively faint at energies above ~ 3 keV). The last column 
in Table 1 lists the average 1.7-3 keV band flux of each observa¬ 
tion, estimated using the best-fit PL model fits to the 3-10 keV 
band data. These values are representative of the continuum flux 
levels in each observation and show that the X-ray continuum 
had the highest and lowest flux levels during Obs. 2 and Obs. 4, 
respectively. Conversely, Obs. 2 (red points in Fig. 1) shows the 
‘weakest’ soft excess, in terms of ratio, while Obs. 4 (blue points 
in the same figure) shows the ‘strongest’ soft excess (i.e. highest 
ratio). 

3.2. The flux-flux plots 

To construct the FFPs, we divide the 0.2 -1.7 keV band into 
11 energy sub-bands, with a width of 0.1 keV from 0.2 up to 
1 ke V. The width increases to 0.2 keV for the next two sub-bands, 
and to 0.3 keV in the final sub-band (1.4 - 1.7 keV), so that the 


signal-to-noise ratio of the resulting light curves in this, and in 
the continuum bands, is larger than 3<x (on average). 

We produced light curves in each energy band (including the 
continuum band) using a bin size, Atbi n , of 1, 4, and 8ks. We 
considered various Afbin values to investigate whether Afbin af¬ 
fects the shape of the resulting FFPs or not. For the given data 
sets, a bin size smaller than 1 ks will result in some data points 
having less than 10 counts in the 1.7-3 keV band. The error bars 
on these data points will be far from Gaussian, and, as a result, it 
would not be possible to fit the FFPs using traditional x 2 statis¬ 
tics. On the other hand, a Afbin larger than 8 ks results in FFPs 
with a few points, and a reduced max/min dynamical range. 

Figure 2 shows the 0.3-0.4, 0.9-1, and 1.4-1.7 vs 1.7-3 keV 
FFPs (top, middle, and bottom rows, respectively) for Obs. 3. 
They are representative of the FFPs of all observations. Left, 
middle, and right columns show the FFPs for Afbin = 1,4, and 
8 ks, respectively. The soft and continuum band count rates are 
highly correlated, as is usual for Seyferts. 

It is customary to fit the FFPs with linear relations in the 
form of: ‘soft band’ flux = ‘hard/reference band’ flux + constant. 
A stable component in the soft band is claimed when the best-fit 
intercept is significantly different from zero. We chose to follow 
a slightly different, and more general, approach in the modelling 
of the FFPs, and we fit them with a power-law (PL) relation in¬ 
stead. If the underlying relation between soft and reference band 
counts is indeed linear, the best-fit model slope should be con¬ 
sistent with unity. If there is a stable component, whose flux is 
smaller than the primary flux, we would expect a flattening of the 
FFP at low reference band fluxes and the detection of a non-zero 
constant, just like when we fit the data with a linear function. 


To investigate the dependence of the FFPs on Af h ,„. we first con¬ 
sidered the ‘high-flux’ data points in these plots, i.e. the points 
with a count rate in the continuum band that is higher than its 
median (the vertical lines in Fig. 2 indicate the median in the 
case of the Obs. 3 plots). In this way, we avoid complications 
that may be associated with a possible flattening at low frequen¬ 
cies, resulting from the presence of a stable spectral component. 

We fitted the high-flux part of the FFPs with a power-law 
(PL) relation of the form 

y = axP , 

where y and x represent the ‘soft’ and the continuum band count 
rates, respectively. The model normalisation, a, and slope, /j, 
were both left to vary freely during the fit. The fits were per¬ 
formed using the Python routine MPFIT 1 , based on the Fortran 
routine LMFIT (More 1978), taking into account the errors on 
the //-axis (i.e. the errors of only the soft energy band counts). 
The solid red lines in Fig. 2 indicate the best-fit PL lines, extrap¬ 
olated to lower fluxes as well. The bottom panels in each plot 
show the best-fit residuals (i.e. the ratio (data-model)/<x). The 
residuals at the high-flux part of the plots are uniformly scat¬ 
tered around zero, and there is no indication of any systematic 
residual trend. 

Figure 3 shows the best-fit PL normalisation and slope values 
for all the FFPs of Obs. 3 (top and bottom panel, respectively). 
Black, red, and blue points indicate the results for the 1, 4, and 
8 ks binned light curves, respectively. The best-fit a and ft values 
of the 4 and 8 ks binned light curves are consistent with each 
other. This is not the case with the 1 ks binned data: the best-fit 


1 https://code.google.com/p/astrolibpy/source/browse/mpfit/mpf 


3.3. The effects of the light curve bin size 
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1.7-3.0 keV Count Rate (Count/s) 1.7-3.0 keV Count Rate (Count/s) 1.7-3.0 keV Count Rate (Count/s) 







1.7-3.0 keV Count Rate (Count/s) 1.7-3.0 keV Count Rate (Count/s) 1.7-3.0 keV Count Rate (Count/s) 

Fig. 2. 0.3-0.4, 0.9 —1, and 1.4-1.7 vs 1.7-3 keV flux-flux plots (top, middle, and bottom rows, respectively) for Obs. 3, obtained with Af bin = 1,4, 
and 8 ks (left, middle, and right columns, respectively). The solid red line indicates the best-fit power-law relation to the ‘high-flux’ data (see text 
for details). The vertical lines indicate the continuum band median. The best-fit residuals are plotted in the lower plan of each plot. Error bars are 
omitted in the residual plots for reasons of clarity . 


a i k. s values are systematically smaller than the respective 4 and 
8 ks values at energies below ~ 1 keV (where the soft excess is 
more pronounced), and the same is true for the best-fit slopes as 

well. 

This can be seen from the plots in the two upper rows in 
Fig. 2. The x— and y —axis range is the same in all panels. It 
is clear that the 4 and 8ks best-fit models (solid red lines in 
the middle and right-hand panels) cannot fit well the high-flux 
part of the 1 ks flux-flux plots (shown in the left panels). We 
observe similar differences between the best-fit ( a\ks,P\ks ) and 
(&4l&ks,/34/$ks) values in the other observations as well. These 
differences strongly suggest that the choice of the time bin size 
does affect the shape of the resulting FFPs. We provide a possi¬ 
ble explanation for the effect of light curve bin size on the FFPs, 
below. 

The 4 and 8 ks FFPs are not a binned version of the 1 ks 
FFPs. The data binning has been performed in the time, and not 
in the flux, domain. A large Afbin value may result in a FFP which 
is distorted and does not correspond to the intrinsic plot, if the 
intrinsic relation between the count rates in two energy bands is 


non linear (i.e. f3 ± 1), and the source is variable. For example, 
supposing that the intrinsic relation between the ‘soft’ and ‘hard’ 
band counts, C s b and Chb, respectively, is in the form of C s b = 
aCy |g, and that both bands are variable. The light curve binning 
implies a binning of the source signal in the time domain in the 
form 

J /V~Ar bm 

Cbin(0 = T—“ C(t)dt , 

A* bin Jt 

where Cbi n stands for ‘binned counts’. If [i mu + 1, then since 

X I+ Afbin / /-V+Afbin 

c ! ^(t)dt + IJ C m (t)dt\ , 

the relation between CHB,bin(0 and CsB,bin(0 will not be identical 
to the intrinsic one. 

IRAS 13224-3809 is a highly variable source in X-rays. 
Gonzalez-Martfn & Vaughan (2012) have estimated its 0.2- 
2 ke V power spectrum. It has a power-law like shape, with a nor¬ 
malization of ~ 10~ 5 and a slope of ~ -2. This implies that the 
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Fig. 3. Best-fit a and /? values (top and bottom panels, respectively), 
obtained by fitting a PL model to the ‘high-flux’ data points (see text 
for details) of the 1, 4, and 8 ks FFPs of Obs. 3 (black filled circles, red 
filled squares, and blue filled triangles, respectively) 


fractional root mean square variability amplitude on timescales 
of 8, 4, and 1 ks is 28%, 20%, and 10%, respectively 2 . So, when 
we bin the light curves using a bin size of 8 and 4 ks, we typ¬ 
ically average data points which scatter around the mean by a 
factor of ~ 0.3 and 0.2, respectively. The average scatter should 
only be ~ 10% in the 1 ks binned light curves. Perhaps then, the 
A/3 differences we observe between the 1 ks and the 4/8 ks flux- 
flux plots are due to the large variability amplitude of the source 
and the intrinsically non-linear flux-flux relations in this source. 
Henceforth, we have decided to study only the 1 ks binned flux- 
flux plots in this work. 


but significant, soft-band variations that are not related to the 
primary source flux. 

The last column in the same table lists the root mean square 
deviation, cr rms , of the data from the model. This is estimated by 


Z 


Uji //model./T 


^model,/ 


1/2 


where N is the number of data points in the FFP, y, and <x,- rep¬ 
resents the observed count rate in the soft energy band and its 
corresponding error, and y mo de\.i represents the model count rate 
in the same band. In practice, cr rms measures the average discrep¬ 
ancy between the best-fit model and the data points, also taking 
the data errors into account (but only of the dependent variable). 
The average data-to-model deviations are ~ 20 - 35% for all 
FFPs of all observations, except Obs. 4, where the deviations 
increase to almost ~ 50% (max). 

Figure 5 shows a plot of the best-fit /ip| c and c values (top 
and bottom panels, respectively) for all observations. The best-fit 
slopes are flatter than 1 at energies below 1 keV, where the soft 
excess is stronger. The only exception is the data from Obs. 2, 
where the PL slopes are close to unity. The best-fit slopes in¬ 
crease towards unity with increasing energy until ~ 0.9 - 1 keV, 
and then are steeper than 1 at higher energies. 

The best-fit c values are consistent with zero in all the FFPs 
of Obs. 5. This is the case with the FFPs of Obs. 4 as well. In 
fact, most of the best-fit c values at energies below 1 keV are 
negative for this observation. On the other hand, the best-fit cs 
are significantly different from zero at energies below ~ 1 keV in 
the FFPs of Obs. 2 and Obs. 3. They are positive in the case of 
the Obs. 1 FFPs as well, although in many cases at a level lower 
than 3<x. 


3.4. The case of a soft-band ‘constant component’ 


3.5. The spectral shape of the soft-band ‘constant 
component' 


The residual’s plots in the left-hand column of Fig. 2 show a data 
‘excess’ above the interpolation of the best-fit PL to low count 
rates. The excess is more pronounced in the top and middle pan¬ 
els (i.e. in the 0.3-0.4 and 0.9-1.0 vs 1.7-3 keV plots), where 
the soft excess is stronger in this source (see Fig. 1). As we have 
argued above, such a flattening in the FFPs at low count rates 
could be due to the presence of a constant soft-band component. 
To quantify its presence, spectral shape, and strength, we fitted 
all the FFPs with a power-law plus constant (PLc) model in the 
form of 

y = «PLc 2 f 8pLc + c. 

We fitted the data using MPFIT. All parameters were left free 
during the fitting. 

The best-fit results are listed in Table 2. The lines within each 
observation list the results for the 0.3-0.4, 0.4-0.5,0.5-0.6, 0.6- 
0.7, 0.7-0.8, 0.9-1, 1-L2, 1.2.-1.4, and 1.4-l-1.7vs 1.7-3 keV 
band FFPs. Figure4 shows the Iks, 0.3-0.4, 0.9-1, and 1.4- 
1.7 vs 1.7-3 keV band flux-flux plots (left, middle, and right 
columns, respectively) for all observations. The dashed lines 
show the best-fit PLc models. The bottom panels in the plots 
show the best-fit residuals, which are uniformly scattered around 
the best-fit models, with no indication of any systematic discrep¬ 
ancies. The PLc model describes the overall trend in the FFPs 
well, and the large x~ reduced values are due to low-amplitude, 

2 By definition, the fractional mean square variability amplitude over a 
given time period, say T, is equal to the integral of the power-spectrum 
density from l/T up to infinity. 


As mentioned above, the non-zero, positive constants detected 
in the FFPs of Obs. 2 and Obs. 3, and perhaps Obs. 1 as well, 
could be indicative of the presence of a separate spectral com¬ 
ponent which is less variable than the primary X-ray source on 
timescales of a few ks. In this case the best-fit c model values 
can be used to study the spectral shape of the constant compo¬ 
nent. For that reason we divided them by the corresponding en¬ 
ergy band width to produce the ’spectrum 1 of this component, in 
photon s^’keV -1 . 

First, we fitted those spectra with an absorbed blackbody 
model (wabs*bbody) by considering only the Galactic absorp¬ 
tion. The best-fit kT\,b turned out to be similar when we fitted 
the three spectra individually, so we repeated the fit by keep¬ 
ing k7bb tied in all of them. The resulting best-fit tempera¬ 
ture was 110 + 1.8 eV, but the overall quality of the fit is poor 
(X 1 /dof = 87.9/29). Then we considered an absorbed power- 
law model (wabs*powerlaw). The best-fit T is steep: ~ 4 for 
Obs. 2 and Obs. 3, and ~ 4.5 for Obs. 1, but the overall quality of 
the fit is worse than before (^ 2 /dof = 137.7/27). 

Then we fitted the spectra with the optxagnf model (Done 
et al. 2012). This is a model for the spectral energy distribution 
of a disc around a rotating SMBH. We first fitted the spectra with 
wabs*optxagnf, assuming Galactic absorption only, and setting 
rcor tied to the same value in all spectra ( rcor , measured in r g , 
sets the radius below which the disc emission emerges as a low- 
temperature, large, optical-depth, Compton-upscattered flux, as 
opposed to a colour-temperature-corrected blackbody at larger 
radii). We ignored the emission of the hard X-ray corona, by 
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Fig. 4. 0.3-0.4, 0.9-1, and 1.4-1.7 vs 1.7-3 keV FFPs (left, middle, and right columns, respectively) for all observations obtained for a bin size of 
1 ks. The dashed red line indicates the best-fit PLc relation. The best-fit residuals are plotted in the lower plan of each plot. The errors are similar 
to those plotted in Fig. 2 but were not plotted for clarity reasons. 
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Table 2. Best-fit PLc models to the FFPs of all observations. 


ffPLc 

A>Lc 

CpLc 

7 

X 

& rms 



Obs. 1 

dof=58 


1.9 ±0.2 

0.69 ± 0.07 

0.07 ± 0.02 

859.8 

0.24 

2.7 ±0.3 

0.81 ±0.06 

0.12 ±0.02 

973.0 

0.25 

2.2 ±0.2 

0.67 ± 0.06 

0.07 ± 0.02 

1028.7 

0.25 

1.7 ±0.2 

0.68 ± 0.07 

0.05 ± 0.02 

880.1 

0.26 

1.5 ±0.2 

0.63 ± 0.07 

0.01 ± 0.02 

591.5 

0.24 

1.3 ±0.2 

0.71 ± 0.07 

0.02 ± 0.01 

602.6 

0.26 

1.5 ±0.2 

0.87 ± 0.07 

0.01 ±0.01 

378.3 

0.27 

1.6 ±0.3 

1.10 ±0.08 

0.01 ± 0.004 

241.7 

0.24 

4.3 ±0.7 

1.56 ±0.08 

0.02 ± 0.002 

331.0 

0.32 

1.1 ±0.2 

1.18 ±0.10 

0.001 ± 0.003 

128.6 

0.20 

1.4 ±0.3 

1.27 ±0.10 

0.003 ± 0.002 

124.4 

0.26 



Obs. 2 

dof=l17 


2.9 ±0.2 

1.02 ±0.04 

0.16 ±0.01 

2003.6 

0.25 

3.1 ±0.2 

1.02 ±0.04 

0.18 ±0.01 

2286.8 

0.26 

3.5 ±0.2 

1.09 ±0.04 

0.19 ±0.01 

2130.3 

0.25 

2.7 ±0.2 

1.06 ±0.05 

0.13 ±0.01 

1553.2 

0.24 

2.6 ±0.2 

1.14 ±0.05 

0.11 ±0.01 

1397.4 

0.25 

1.8 ±0.1 

0.99 ± 0.05 

0.06 ± 0.01 

922.7 

0.22 

1.9 ±0.2 

1.17 ±0.06 

0.05 ± 0.01 

669.0 

0.21 

1.5 ±0.1 

1.22 ±0.06 

0.03 ± 0.004 

453.7 

0.21 

2.0 ±0.2 

1.27 ±0.05 

0.03 ± 0.004 

471.6 

0.21 

1.1 ±0.1 

1.19 ±0.06 

0.01 ± 0.003 

293.1 

0.20 

1.2 ±0.1 

1.28 ±0.07 

0.01 ± 0.003 

233.6 

0.14 



Obs. 3 

dof=l17 


2.3 ±0.2 

0.80 ± 0.04 

0.08 ± 0.01 

2806.7 

0.35 

2.2 ±0.1 

0.71 ±0.04 

0.06 ± 0.01 

3437.1 

0.35 

2.8 ±0.2 

0.85 ± 0.04 

0.10 ±0.01 

3283.3 

0.35 

2.3 ±0.2 

0.87 ± 0.04 

0.06 ± 0.01 

2491.0 

0.35 

2.0 ±0.2 

0.90 ± 0.04 

0.05 ± 0.01 

1822.8 

0.33 

1.8 ±0.2 

0.96 ± 0.05 

0.04 ± 0.01 

1241.2 

0.30 

1.1 ±0.1 

0.82 ± 0.05 

0.01 ±0.01 

922.1 

0.30 

1.4 ±0.2 

1.12 ±0.06 

0.02 ± 0.003 

550.3 

0.26 

2.1 ±0.2 

1.21 ±0.05 

0.01 ± 0.003 

545.0 

0.28 

1.1 ±0.2 

1.24 ±0.07 

0.01 ± 0.002 

347.1 

0.28 

1.1 ±0.2 

1.21 ±0.07 

0.004 ± 0.002 

268.6 

0.24 



Obs. 4 

dof=107 


3.4 ±0.2 

0.77 ± 0.03 

-0.01 ±0.01 

2815.9 

0.58 

3.3 ±0.2 

0.72 ± 0.03 

-0.03 ± 0.01 

3225.0 

0.58 

3.2 ±0.2 

0.73 ± 0.03 

-0.02 ± 0.01 

3015.6 

0.55 

2.7 ±0.2 

0.78 ± 0.03 

-0.01 ±0.01 

7119.4 

0.56 

2.4 ±0.2 

0.80 ± 0.03 

-0.01 ± 0.005 

1804.5 

0.56 

2.3 ±0.2 

0.90 ± 0.04 

-0.001 ± 0.004 

1275.4 

0.56 

1.6 ±0.2 

0.88 ± 0.04 

-0.01 ± 0.003 

998.8 

0.72 

1.3 ±0.2 

1.00 ±0.05 

-0.002 ± 0.001 

575.9 

0.71 

2.0 ±0.4 

1.23 ±0.07 

0.004 ± 0.001 

301.7 

0.42 

0.9 ±0.2 

1.17 ±0.09 

2.1E-4± 0.001 

245.9 

0.51 

0.9 ±0.2 

1.14 ±0.08 

2.7E-5 ± 0.001 

198.6 

0.54 



Obs. 5 

dof=l17 


3.2 ±0.1 

0.79 ± 0.02 

0.004 ± 0.01 

1909.0 

0.31 

3.6 ±0.1 

0.80 ± 0.02 

-0.001 ±0.01 

2173.3 

0.32 

3.8 ±0.1 

0.83 ± 0.02 

0.01 ±0.01 

2203.7 

0.32 

2.7 ±0.1 

0.81 ±0.02 

0.005 ± 0.01 

1553.7 

0.30 

2.6 ±0.1 

0.89 ± 0.02 

0.02 ± 0.01 

1291.3 

0.29 

2.2 ±0.1 

0.93 ± 0.02 

0.01 ± 0.004 

1070.9 

0.33 

1.9 ±0.1 

1.00 ±0.03 

0.01 ± 0.003 

688.4 

0.30 

1.5 ±0.1 

1.09 ±0.03 

0.01 ± 0.002 

441.8 

0.29 

2.3 ±0.1 

1.27 ±0.03 

0.01 ± 0.002 

389.6 

0.24 

1.4 ±0.1 

1.32 ±0.04 

0.004 ± 0.001 

279.8 

0.27 

1.0 ±0.1 

1.16 ±0.04 

0.001 ± 0.001 

205.2 

0.23 



Energy (keV) 

Fig. 5. Best-fit parameters /? PLc (top) and c (bottom), obtained by fitting 
a PLc model to the 1-ksec binned light curves for Obs. 1 (black filled 
squares), Obs. 2 (black open squares), Obs. 3 (red filled circles), Obs. 4 
(red open circles), and Obs. 5 (blue filled triangles).The model parame¬ 
ters, obtained by fitting a PLc relation to simulated FFPs, are also plot¬ 
ted for the high- and low-flux case models we considered (magenta- 
filled diamonds and magenta-filled hexagons, respectively). The model 
parameters, in the case of a variable primary in flux and spectral slope 
without adding an excess component, are also plotted for the high- 
flux and low-flux cases (open magenta diamonds and hexagons, respec¬ 
tively; see Sec. 3.6 for details). 

keeping the model parameter^/?/ frozen to 10 6 , since we did 
not take the full band spectrum of the source into considera¬ 
tion. To constrain the fit as much as possible, we kept the black 
hole (BH) mass frozen to 10 7 M e (e.g. Zhou & Wang 2005; Em- 
manoulopoulos et al. 2014). We also considered a maximally 
rotating BH by freezing the spin parameter, astar, to its maxi¬ 
mum value of 0.998 in all spectra. The best-fit Eddington ratio 
(log(L/LEdd)) turned out to be similar for Obs. 1 and 3, and for 
that reason we repeated the fit by keeping it tied for these two 
observations. We also kept kT e and r (the electron temperature 
and optical depth for the soft Comptonization component) tied 
to the same values in all spectra. 

The quality of the fit improves significantly (x 1 1 dof = 
52.6/28), but this is still not a statistically accepted fit. We ob¬ 
tained a 3<x lower limit on the coronal radius of rcor > 1.93 r g . 
The best-fit Eddington ratio is log(L/LEdd) = -1.36'[| | s | for 
Obs.l and 3, and log(L/LEdd) = -1 .lO'Jj/Z for Obs. 2. Although 
the best-fit required different accretion ratios for the Obs. 1 and 
3 and Obs. 2 spectra, the resulting best-fit values are consistent 
with each other within the errors. The best-fit electron tempera¬ 
ture and optical depth were kT e = 0.15*o[^keV and r = 22^'°, 
respectively. 

3.6. The slope of the flux-flux plots 

The PLc best fit slope values of the FFPs at energies below 
~ 1 keV are less than unity in most observations. This result 
suggests the presence of intrinsic spectral variations in the con¬ 
tinuum emission of the source. To investigate this issue further, 
we used the SAS tool efluxer and estimated the source flux 
(in erg/cm 2 /s) in all the energy bands of all observations. We 
also fitted the Obs. 1-5 spectra in the 1.7-3 keV band with an ab¬ 
sorbed PL model. We computed the best-fit PL flux in the same 
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Fig. 6. F ex C ess,obs plotted as a function of FpL.obs, for Obs.l, 2, 3, 4, and 5 
(black, red, cyan, blue, and magenta squares, respectively). The errors 
on the points were calculated using the efluxer flux error and the error 
on the best-fit slope, and normalisation of the PL model fit to the 1.7-3 
keV data. The dashed lines indicate the best-fit PL model fits to the data. 


bands and we calculated the difference between the observed and 
the PL model flux. This should be representative of the ‘intrinsic’ 
soft-band flux alone, F exc ess,obs- Figure 6 shows F excesSj obs versus 
the 1.7-3 keV flux (which should be representative of the X- 
ray primary flux, Fpr.obs- for the 0.2-0.3, 0.5-0.6, 0.8-0.9, and 
1.2-1.4 bands. (These plots are representative of the plots in all 
energy bands.) . The observed excess flux in the 1.2-1.4 keV 
band (and in general at energies above 1 keV) is almost an or¬ 
der of magnitude smaller than the one at energies below ~ 1 
keV, and it is hard to constrain. In fact, it may even be zero (in 
agreement with Fig. 1, which shows that the soft excess flux at 
energies above ~ 1 keV is minimal). 

We fitted the data in all energy bands with a PL model, with 
the slope fixed to 0.46. This is the best-fit slope value we found 
when we fitted the data in panel (a) of Fig. 6 in Chiang et al. 
(2015). These data show the relation between the 0.1-100 keV 
flux of the reflection component (which may be responsible for 
the soft excess) and the primary flux (in the same band) in IRAS 
13224-3809. This PL model fits the “Ferness,obs vs FpL.obs” data 
in all bands reasonably well. 

We then produced synthetic energy spectra (using the XSPEC 
command fakeit), assuming a PL model and taking the Galac¬ 
tic absorption into account, in which the spectral slope of the 
continuum varies with the continuum flux according to F = 
5.73(Vp L '. This is the best-fit model to the slope and normalisa¬ 
tion of the PL fits to the 1.7-3 keV band data of each observa¬ 
tion, and it it is in good agreement with the spectral variability 
law (T oc Apj 08 ) found by Sobolewska & Papadakis (2009) .We 
considered two cases for the model variability: (a) We first tried 
to reproduce the source behaviour during Obs. 2, which shows 
the highest flux, which meant we considered a mean continuum 
flux equal to Fpr.obst and an Npr variability of a factor of 30, 
similar to the observed one (the high-flux case hereafter), (b) We 
did the same, but tried to model the variability seen in Obs. 4, 
which shows the lowest flux and an Npp variation flux variation 
of ~ 50 (the low-flux case hereafter). The high-flux case implies 
T variation between 1.95 and 2.95, while the spectral slope varies 
between 1.6 and 2.55 in the low-flux case. 


We measured the flux of the simulated spectra in the soft- 
bands we consider in this work, as well as in the 1.7-3 keV 
band (FpL.simui)- We then added an ‘excess’ flux to the soft-band 
fluxes, which we calculated using Fpi,.si mu i and the best-fit PL 
models to the observed ‘F exC ess,obs ~ Fpi ohs ’ plots in each band. 
Next, we plotted the total soft-band flux, F SO ft iS i mu i, as a func¬ 
tion of FpL.simui for all the Npp values. The resulting ‘synthetic’ 
flux-flux plots are plotted in Fig. 7. We fitted the simulated FFPs 
with a PLc model, exactly as we did with the observed FFPs. 
The dashed and solid lines in Fig. 7 show these best-fit models 
(for the high- and low-flux cases, respectively). Interestingly, we 
found that the constant c in the PLc model fits was not equal to 
zero. In other words, a flattening in the FFPs is expected, even 
when there is no constant soft component, but the soft-excess 
component is variable, albeit in a fashion correlated with the PL 
continuum. 

The best-fit ySpLc.mod and cpL C ,mod model values are plotted 
in Fig. 5. To transform the best-fit cpL C ,mod values (currently in 
erg/cm 2 /s units) into count/s, we used the HEASOFT web tool 
PIMMS, assuming that the broad-band soft-excess component in 
IRAS 13224-3809 is approximated by a blackbody model, with 
kT = 0.1 keV. (Such a blackbody component does approximate 
the soft excess in this object in very broad terms, as we verified 
by fitting a PL plus bbody model to the data shown in Fig. 1.) 

The high-flux case and low-flux case /ippc.mod predictions 
bracket the observed values at energies below ~ 0.9 keV. At 
higher energies, both case models predict FFPs with slopes flat¬ 
ter than the observed ones. This is probably because the soft- 
excess contribution to the observed flux is minimal at high en¬ 
ergies. For example, the empty diamonds and hexagons in the 
top panel of Fig. 5 indicate the /?pLc,mod values when we only 
consider the PL continuum with variable T, i.e. without adding 
Fexcess, for both cases. The agreement with the observed slopes 
is much better in this case. 

The model constants cpLc.mod are smaller than the observed 
constants, both in the high- and low-flux cases. This indicates 
that, although a flattening at low fluxes is expected in the case of 
a primary variable in norm and shape and of a variable soft-band 
component, a stable component may still exist. The spectrum 
of this component can be fitted well with an optxagnf model 
(X 2 /dof = 12.7/10) by fixing all the parameters to the best- 
fit values reported in Sec. 3.5 for Obs. 2 and only leaving the 
log(L/LEdd) as a free parameter. The new best-fit value for the 
Eddington ratio is ~ 0.05, which is - 2 times less than the pre¬ 
vious value for this observation. 


4. Discussion and conclusions 

We have presented a timing analysis of the five XMM-Newton 
archival observations of the NLS1 Galaxy IRAS 13224 - 3809. 
Our main aim was to study its X-ray variability at energies below 
~ 1.5 keV. 

We produced FFPs using the data in 11 energy bands be¬ 
tween 0.2 and 1.7 keV and in the 1.7-3 keV band data, which we 
assumed is representative of the X-ray primary ‘continuum’. We 
considered various bin sizes for the light curves, and we found 
that the bin size choice is important because it affects the shape 
of the resulting plots. This is probably the case in sources that 
show fast variations on short timescales, and where the intrin¬ 
sic FFPs are non-linear. In these cases, the choice of the shortest 
possible bin size is preferred. 

Our main results, using the 1-ks binned light curves, are sum¬ 
marised below: 
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Fig. 7. Plots of F so f tiS i mu i as a function of FpL.simui in the high-flux and 
low-flux variability cases (red circles and black squares, respectively). 
The FPPs were fitted with a PLc model (dashed red line and solid black 
line for the high- and low-flux cases, respectively). 


- We fitted the FFPs with a PLc model. Positive constants are 
detected at energies below 1 keV (where the soft excess is 
more pronounced in this source) in Obs. 1, 2, and 3. All the 
best-fit values are consistent with zero in the case of Obs. 5. 
We also detected negative constants in Obs. 4, which are all 
negative at all energies less than 1 keV, although they are 
consistent with zero (at 3<x). 

- The best-fit slopes are significantly flatter than 1 at energies 
below ~ 1 keV, except in Obs. 2, where they are slightly 
steeper than unity. This indicates that the intrinsic flux-flux 
relation is not linear and points to intrinsic spectral variabil¬ 
ity. 

Strictly speaking, the PLc model is not statistically accepted. 
Nevertheless, the model represents the general trends in the flux- 
flux plots well. The residual plots do not show any large-scale, 
systematics residuals, which is indicative of the presence of an 
extra, broad-band model component. The high best-fit x 1 values 
are due to random data fluctuations around the best-fit models, 
which have an amplitude of ~ 20 - 55% (of the model value) 
in all observations. The observed range of variations in all ob¬ 
servations is comparable to 5-10 in the soft energy bands. The 
PLc model, therefore, does take account of most of the observed 
variations, and there is just a scatter, of ~ 0.2 - 0.5, that is left in 
the residual’s plots. This implies that short-amplitude, fast varia¬ 
tions in the soft energy bands exist and do not depend on the hard 
band flux. However, the study of the causes of these variations is 
beyond the scope of the present work. 

The flattening at low-flux rates that we detect in Obs. 1, 2, 
and 3 at energies below ~ 0.9-1 keV could imply the presence of 
a separate component that is not variable on timescales of a few 
ks. However, the resulting spectrum of this component cannot be 
fitted by either a blackbody, power-law, or by optxagnf (which is 
a more realistic model for the accretion disc emission in AGN ( 
Done et al. 2012)). 

For almost all observations, the best-fit PLc slopes are not 
equal to unity at all energies. This implies significant intrinsic 
spectral variations. We investigated the case of a primary con¬ 
tinuum, which is variable in flux, and spectral slope as F oc A' 44 
(in agreement with Sobolewska & Papadakis 2009) as well as 


a soft band excess flux, which is variable and which positively 
correlated with the continuum according to the relation F excesfi oc 
A 0 46 . We find that the slope of the resulting model FFPs de- 

primary r ° 

pends on energy in a similar way to what we observed. The slope 
of the model FFPs at energies above ~ 0.9 keV is similar to the 
slope of the observed plots, but only if the soft-band component 
does not contribute to these energies. In addition, the model FFPs 
flatten at low fluxes, resulting in the detection of a positive con¬ 
stant when fitted by a PLc model, although no constant soft-band 
component is present. 

The constant values that we obtained from the simulated 
FFPs are lower than the ones we observed. This may indicate 
the presence of a separate, constant component. The 0.2-1 keV 
band flux of this component should be ~ 15% of the total ob¬ 
served flux during Obs. 2. In the case of Obs 2, however, this 
component is broadly consistent with the energy spectrum emit¬ 
ted by an accretion disc around a maximally rotating BH of mass 
~ 10 7 M© with an accretion rate of ~ 0.05 of the Eddington limit. 
This is a rather low value for NLS1 galaxies, which are believed 
to accrete at much higher rates. At radii r < r cm ona = 2r g (at 
least), the energy is dissipated by Compton up-scattering of seed 
photons at the disc temperature at r - r corona , off electrons with 
a temperature of ~ 0.15 keV and optical depth of ~ 22. How¬ 
ever, the extrapolation of best-fit optxagnf model to the ultravi¬ 
olet band at 2310 A gives a flux of 1.2 X 10~ 16 erg cm~ 2 s _1 A _1 
in the case of Obs. 2. This is just ~ 4% of the observed OM 
UVM2 flux, which is 3.07 x 10 -15 ergcm _2 s _1 A _1 , as reported 
by Chiang et al. (2015). 

For this reason, we believe that the positive constants we de¬ 
tect in the FFPs are mainly due to the intrinsic spectral vari¬ 
ations. Despite the agreement between the model and the ob¬ 
served FFPs’ slope, it is difficult to explain the non-detection 
of positive constants in the Obs. 5 FFPs. Even if all the posi¬ 
tive constants are caused by intrinsic, complicated spectral vari¬ 
ations, their non-detection in Obs. 5 implies that the source op¬ 
erates in a different way in this case (and yet the flux-flux slopes 
are similar to those observed in the other observations). It is even 
more difficult to explain the negative constants that we observe 
in the Obs. 4 FFPs. One possibility is that the source is affected, 
at least occasionally, by extra, intrinsic absorption. This cannot 
be neutral; an extra neutral absorber can result in ~ zero, but 
not negative, constants. In addition, the slopes will be signif¬ 
icantly flatter than the observed ones in the low energy FFPs. 
On the other hand, an ionized absorber that is flux-dependent, in 
such a way that the absorption is stronger at low fluxes, may be 
able to explain the presence of negative constants in the FFPs. 
In practice, a warm absorber can vary in column density, cover¬ 
ing fraction, and/or ionization parameter. The study of this kind 
of absorber is beyond the scope of the present work, since the 
modelling of its effects is highly complicated and not well con¬ 
strained if all the parameters of the absorber are variable. 

In summary, our results support the hypothesis that most of 
the soft excess flux at energies below ~ 0.9 keV is due to X-ray 
reflection in IRAS 13224-3809. This agrees with the results of 
Chiang et al. (2015) and with those of Emmanoulopoulos et al. 
(2014) and Kara et al. (2013), who interpret the observed soft 
band time-lags in the same scenario. The soft excess is respond¬ 
ing to the primary X-ray variations, although with a smaller am¬ 
plitude (as expected for a smeared component). At the same 
time, the primary slope steepens with increasing flux (which is 
again agrees with the results of Chiang et al. 2015). We cannot 
exclude the presence of an extra, stable soft-band component, 
which can be well fitted by a disc plus a warm Comptonization 
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medium emission. Its contribution to the observed flux, and at 
energies below 1 keV, should be less than ~ 15%. 

The study of the FFPs can yield interesting results, not only 
when positive constants are detected (which may not correspond 
to a constant component), but also when none or even negative 
constants are detected. This is the case with two observations 
of IRAS 13224-3809, and they may indicate the presence of an 
intrinsic warm and variable absorber in the source. 
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